# 1.1 Read region labels (World Bank country & lending groups)
lbl_data <- read_xlsx(here("lbl_regions.xlsx"))
# ‣ Columns: Entity (OWID name), Code (ISO-3), Region (WB group)
# 1.2 Download per-capita electricity mix
grid_data <- read_csv(
"https://ourworldindata.org/grapher/per-capita-electricity-fossil-nuclear-renewables.csv?v=1&csvType=full&useColumnShortNames=true"
) |>
rename(
fossil_pc = per_capita_fossil_generation__kwh_chart_per_capita_electricity_fossil_nuclear_renewables,
nuclear_pc = per_capita_nuclear_generation__kwh_chart_per_capita_electricity_fossil_nuclear_renewables,
renewable_pc = per_capita_renewable_generation__kwh_chart_per_capita_electricity_fossil_nuclear_renewables
)
# 1.3 Download carbon intensity (g CO₂e/kWh)
carbon_intensity <- read_csv(
"https://ourworldindata.org/grapher/carbon-intensity-electricity.csv?v=1&csvType=full&useColumnShortNames=true"
)
# 1.4 Merge, fallback, filter & compute generation shares
grid_data_ci <- grid_data |>
# join in Region; if missing, use OWID name as region
left_join(lbl_data |> select(-Entity), by = "Code") |>
mutate(Region = if_else(is.na(Region), Entity, Region)) |>
# attach CO₂ intensity and drop missing
left_join(
carbon_intensity |>
select(Entity, Year, co2_intensity = co2_intensity__gco2_kwh),
by = c("Entity", "Year")
) |>
drop_na(co2_intensity) |>
# pick most recent per Entity
group_by(Entity) |>
arrange(desc(Year), .by_group = TRUE) |>
slice(1) |>
ungroup() |>
# compute percent shares of generation mix
mutate(
total_generation = fossil_pc + nuclear_pc + renewable_pc,
fossil_generation = 100 * fossil_pc / total_generation,
nuclear_generation = 100 * nuclear_pc / total_generation,
renewable_generation = 100 * renewable_pc / total_generation
) |>
# re-attach co2_intensity (to preserve Year) and fix Europe code
select(-co2_intensity) |>
left_join(
carbon_intensity |>
select(Entity, Year, co2_intensity = co2_intensity__gco2_kwh),
by = c("Entity", "Year")
) |>
mutate(Code = if_else(Entity == "Europe", "EU", Code))Compute Oxygen Emissions
Part 1: Data Acquisition & Preparation
What we do: 1. Load a hand‐coded World Bank Country→Region mapping. 2. Download per-capita electricity generation mixes from OWID. 3. Download grid carbon intensity (g CO₂e /kWh) from OWID. 4. Merge everything and keep only the most recent year with a non-missing CI per country.
Why: - We need country-level grid carbon intensities and generation shares to scale our oxygen LCA. - Regions are used later for delivery‐mix weighting.
Part 2: Define Emissions for Each Delivery System
Goal: Break down system‐level CO₂e into legs (“prod”, “transp”, “otro”), scale to 100% purity, and convert to g CO₂e per pure liter.
Key logic: - We measure a total GWP per liter at some purity → allocate by % to each leg → “undo” purity → convert to grams.
# Helper fn: compute component‐level GWP at 100% purity
make_system_emissions <- function(
components, # leg names
percent_emissions, # numeric vector summing to 100
total_kg_GWP_1L, # observed kg CO2e / L at measured purity
purity_factor # fraction (e.g. 0.995)
) {
tibble(
component = components,
percent_emissions = percent_emissions
) |>
mutate(
# raw kg per component at measured purity
kg_GWP_1L_raw = percent_emissions / 100 * total_kg_GWP_1L,
# scale to 100% purity
kg_GWP_1L = kg_GWP_1L_raw / purity_factor,
# grams per pure L
g_GWP_1L = kg_GWP_1L * 1000
) |>
select(component, percent_emissions, kg_GWP_1L, g_GWP_1L)
}
# System 1: Bulk LOX
system1_lox <- make_system_emissions(
components = c("total", "prod", "transp", "otro"),
percent_emissions = c(100, 82, 18, 0),
total_kg_GWP_1L = 0.0001700,
purity_factor = 0.995
)
# System 2: Cylinder LOX
system2_cylox <- make_system_emissions(
components = c("total", "prod", "transp", "otro"),
percent_emissions = c(100, 34, 36, 30),
total_kg_GWP_1L = 0.0004110,
purity_factor = 0.995
)
# System 3: PSA_local
system3_psal <- make_system_emissions(
components = c("total", "prod", "transp", "otro"),
percent_emissions = c(100, 99, 0, 1),
total_kg_GWP_1L = 0.0000844,
purity_factor = 0.93
)
# System 4: Oxygen concentrator
system4_oxconc <- make_system_emissions(
components = c("total", "prod", "transp", "otro"),
percent_emissions = c(100, 99.5, 0, 0.5),
total_kg_GWP_1L = 0.000123,
purity_factor = 0.93
)
# System 5: Composite PSA cylinder
# - prod & otro from PSA local (93%)
# - transp from cylinder LOX (undo 99.5%→93%)
purity_cylox <- 0.995
purity_psal <- 0.93
psal_base <- system3_psal |>
filter(component %in% c("prod", "otro")) |>
select(component, kg_GWP_1L, g_GWP_1L)
cylox_transp_adj <- system2_cylox |>
filter(component == "transp") |>
mutate(
raw_kg = kg_GWP_1L * purity_cylox,
kg_GWP_1L = raw_kg / purity_psal,
g_GWP_1L = kg_GWP_1L * 1000
) |>
select(component, kg_GWP_1L, g_GWP_1L)
components_only <- bind_rows(psal_base, cylox_transp_adj)
new_total_kg <- sum(components_only$kg_GWP_1L)
new_total_g <- sum(components_only$g_GWP_1L)
system5_psac <- bind_rows(
tibble(
component = "total",
percent_emissions = 100,
kg_GWP_1L = new_total_kg,
g_GWP_1L = new_total_g
),
components_only |>
mutate(percent_emissions = 100 * kg_GWP_1L / new_total_kg)
) |>
arrange(factor(component, levels = c("total", "prod", "transp", "otro")))Liquid Oxygen
Cylinder of LOX
PSA Piped
Oxygen Concentrator
PSA Cylinder
Part 3: Calculate Emissions per Country per Delivery System
Objective: For each country and each of our five delivery systems, compute the total CO₂e (g) per pure-L of O₂ by: 1. Combining the per-system component tables into one long format. 2. Cross-joining with each country’s grid carbon intensity (CI) and scaling the “production” leg proportionally to the country’s CI. 3. Re-normalizing the component shares to sum to 100%, summing the grams CO₂e across legs, pivoting to a wide layout, and attaching these values back to the country CI table.
# 3.1 Combine all five system breakdowns into one long table
# - 'systems' is a named list of tibbles: one per delivery system
# - imap_dfr() binds them row-wise and adds a 'system' column
systems <- list(
system1_lox = system1_lox,
system2_cylox = system2_cylox,
system3_psal = system3_psal,
system4_oxconc = system4_oxconc,
system5_psac = system5_psac
)
systems_long <- imap_dfr(
systems,
~ mutate(.x, system = .y) # .x = each tibble; .y = its name
)
# 3.2 Cross-join each country with every system component,
# then scale only the “prod” (production) leg according to the ratio:
# country_CI / ONTARIO_BASELINE_CI
ONTARIO_BASELINE_CI <- 28 # g CO₂e per kWh from Ontario LCA
emissions_components <- grid_data_ci %>%
select(Entity, co2_intensity) %>% # keep only country name & its CI
crossing(systems_long) %>% # replicate each component for each country
mutate(
# if this row is the production leg, scale g_GWP_1L by CI ratio
g_GWP_1L = if_else(
component == "prod",
g_GWP_1L * (co2_intensity / ONTARIO_BASELINE_CI),
g_GWP_1L
)
) %>%
group_by(Entity, system) %>%
mutate(
# re-compute each leg’s % share so that prod + transp + otro = 100
percent_emissions = 100 * g_GWP_1L / sum(g_GWP_1L[component != "total"])
) %>%
ungroup()
# 3.3 Summarize to one total g CO₂e per pure-L for each country × system,
# then pivot to wide format and join back to the country CI table.
emissions_totals <- emissions_components %>%
filter(component != "total") %>% # drop the original 'total' row
group_by(Entity, system) %>%
summarise(
total_g_GWP_1L = sum(g_GWP_1L), # grams CO₂e per pure-L
.groups = "drop"
)
emissions_wide <- emissions_totals %>%
# prefix system names with "GWP_" for clarity in column names
mutate(system = paste0("GWP_", system)) %>%
pivot_wider(
names_from = system,
values_from = total_g_GWP_1L
)
# Attach the per-system GWP columns to the original country CI table
country_emissions <- grid_data_ci %>%
left_join(emissions_wide, by = "Entity")Part 4: Calculate emission 1 L/min per country
Objective Blend the five delivery‐system GWP factors with region‐specific usage shares to derive a single g CO₂e per pure-L of O₂ per minute for every country.
Key Assumptions
- Regional Shares The
oxygen_distributiontable reflects typical delivery-system splits by region (Supplemental Table X). Shares sum to 100% in each region. - Continental Aggregation Rows labeled “Africa”, “Europe”, etc., are not raw data but aggregated averages of their constituent sub-regions.
Steps
- Define regional delivery-system mix shares (Supplemental Table X).
- Compute weighted average g CO₂e per pure-L using those shares.
- Compute continent-level means for “whole” regions and override.
- Save
oxemissions_bycountry.csv.
4.1 Delivery-mix tribble
# 4.2 Prepare the regional mix shares for weighting ---------------------------
# Append “_share” to each system column so it lines up with the GWP columns
oxygen_mix <- oxygen_distribution %>%
rename_with(~ paste0(.x, "_share"), -Entity)
# List of the per-system GWP columns produced in Part 3
system_cols <- c(
"GWP_system1_lox",
"GWP_system2_cylox",
"GWP_system3_psal",
"GWP_system4_oxconc",
"GWP_system5_psac"
)
# Corresponding list of the share columns from oxygen_mix
share_cols <- c(
"system1_lox_share",
"system2_cylox_share",
"system3_psal_share",
"system4_oxconc_share",
"system5_psac_share"
)
# 4.3 Compute weighted average emissions per pure-L for each country -----------
country_1l <- country_emissions %>%
# bring in the regional system shares
left_join(oxygen_mix, by = c("Region" = "Entity")) %>%
# rowwise allows us to compute a per-row sum across multiple columns
rowwise() %>%
mutate(
# sum(GWP_i * share_i) for i = 1…5, then divide by 100 because shares are %
g_CO2e_1L = sum(
c_across(all_of(system_cols)) * c_across(all_of(share_cols)),
na.rm = TRUE
) /
100
) %>%
ungroup()
# 4.4 Compute continent-level averages and override continent rows ------------
# Which Entity names represent whole continents
continent_names <- c("Africa", "Asia", "Europe", "Oceania", "South America")
# 1) Exclude continent rows, map each sub-region to its continent,
# then compute the mean g_CO2e_1L for each continent
continent_avgs <- country_1l %>%
filter(!Entity %in% continent_names) %>%
mutate(
continent = case_when(
Region %in% c("SSA", "MENA") ~ "Africa",
str_detect(Region, regex("Asia", ignore_case = TRUE)) ~ "Asia",
str_detect(Region, regex("Europe", ignore_case = TRUE)) ~ "Europe",
str_detect(Region, regex("Oceania", ignore_case = TRUE)) ~ "Oceania",
str_detect(Region, regex("Latin America", ignore_case = TRUE)) ~
"South America",
TRUE ~ NA_character_
)
) %>%
filter(!is.na(continent)) %>%
group_by(continent) %>%
summarise(
g_CO2e_1L = mean(g_CO2e_1L, na.rm = TRUE),
.groups = "drop"
)
# 2) Join those continent means back to the full table
# and overwrite the continent rows with their average
country_final <- country_1l %>%
left_join(continent_avgs, by = c("Entity" = "continent")) %>%
mutate(
g_CO2e_1L = if_else(
Entity %in% continent_names,
g_CO2e_1L.y, # continent mean
g_CO2e_1L.x # original country value
)
) %>%
# drop the helper columns and any zero‐emission rows
select(-g_CO2e_1L.x, -g_CO2e_1L.y) %>%
filter(g_CO2e_1L > 0)
# 4.5 Render the final table as an interactive DataTable -----------------------
DT::datatable(
country_final,
extensions = 'Buttons',
options = list(
pageLength = 25, # show 25 rows at a time
scrollX = TRUE, # horizontal scrolling
dom = 'Bfrtip', # show Buttons, filter, table, info, pagination
buttons = c('copy', 'csv', 'excel', 'pdf')
),
caption = '### Baseline O₂ Emissions per L by Country'
)Part 5: Projected Grid Carbon Intensity Scenarios for 2030
In this section we augment our baseline country data with three hypothetical future grid carbon‐intensity metrics (g CO₂e /kWh), under differing decarbonization pathways:
- 2030 Renewable Target Scenario (
co2_intensity_2030)- Assumes each country achieves its stated 2030 renewable energy share
- Nuclear share remains fixed at today’s level
- Fossil reduction is entirely backfilled by renewables
- Half-Fossil Scenario (
co2_intensity_half)- Assumes current fossil generation share is cut in half
- Nuclear share remains constant
- Displaced fossil generation replaced by renewables
- Clean-Grid Scenario (
co2_intensity_clean)- Assumes 100 % elimination of fossil generation
- Nuclear share remains at baseline
- Residual energy supplied solely by renewables
Key Assumptions & Parameters
- Nuclear emission factor: 12 g CO₂e /kWh (IPCC AR5)
- Pure-fossil & pure-renewable factors: Median carbon intensities of grids with > 99 % fossil or > 99 % renewables
- Imputation Rules for Missing Targets:
- European countries with no target → use “Europe” (EU) average
- South America → inherits mean of all Latin America sub-regions
- Remaining Latin America gaps → filled with the South America value
We compute these three new intensity columns in our target_intensities table, then feed them into Part 6’s per-liter emissions calculations.```
# 5.1 Read & merge 2030 RES targets
targets_renwbl <- read_xlsx(here("targets_download.xlsx"), sheet=4) |>
select(Code=country_code, target_perc_2030=res_share_target) |>
mutate(Code=if_else(Code=="XKX","OWID_KOS",Code))
target_2030 <- country_final |>
left_join(targets_renwbl, by="Code") |>
mutate(
target_perc_2030 = if_else(
is.na(target_perc_2030) & str_starts(Region,"Europe"),
first(target_perc_2030[Entity=="Europe"], NA_real_),
target_perc_2030
),
target_perc_2030 = if_else(
Entity=="South America",
mean(target_perc_2030[str_starts(Region,"Latin America")], na.rm=TRUE),
target_perc_2030
),
target_perc_2030 = if_else(
is.na(target_perc_2030) & str_starts(Region,"Latin America"),
first(target_perc_2030[Entity=="South America"], NA_real_),
target_perc_2030
)
) |>
filter(!is.na(Code))
# 5.2 Define pure emission factors
ef_nuclear <- 12
ef_fossil <- grid_data_ci |> filter(fossil_generation>99) |> summarise(median(co2_intensity,na.rm=TRUE)) |> pull()
ef_renew <- grid_data_ci |> filter(renewable_generation>99)|> summarise(median(co2_intensity,na.rm=TRUE)) |> pull()
# 5.3 Compute scenario intensities
target_intensities <- target_2030 |>
mutate(
co2_intensity_2030 = if_else(
is.na(target_perc_2030), NA_real_,
{
rg <- pmin(100 - nuclear_generation, pmax(renewable_generation, target_perc_2030))
fg <- 100 - nuclear_generation - rg
fg/100*ef_fossil + rg/100*ef_renew + nuclear_generation/100*ef_nuclear
}
),
co2_intensity_half = {
fg <- fossil_generation * 0.5
rg <- 100 - nuclear_generation - fg
fg/100*ef_fossil + rg/100*ef_renew + nuclear_generation/100*ef_nuclear
},
co2_intensity_clean = {
rg <- 100 - nuclear_generation
rg/100*ef_renew + nuclear_generation/100*ef_nuclear
}
)Part 6: Apply Scenario Intensities & Compile Final Emissions Table
In this final stage, we take the three grid-intensity projections from Part 5 and feed them through our delivery-system LCA pipeline (as in Part 3), apply the regional mixes from Part 4, and then stitch everything together into one comprehensive dataset.
Workflow Overview: 1. Inject Scenario CI - For each scenario (2030 target, half-fossil, clean-grid), replace the baseline co2_intensity in grid_data_ci with the scenario-specific column (co2_intensity_2030, co2_intensity_half, or co2_intensity_clean). 2. Recompute System-Level Emissions - Cross-join with systems_long to get every country × delivery-system component. - Scale the production leg by (scenario CI / Ontario baseline CI). - Recalculate component shares so that prod + transp + other = 100 %. 3. Aggregate to Country × System - Sum the grams CO₂e per pure-L for each system, then pivot to wide format (e.g. GWP_system1_lox, …). 4. Compute 1 L/min Emissions per Country - Apply the regional delivery-mix shares to each system’s wide table to get a weighted g_CO2e_1L_<scenario> for every country. 5. Merge Baseline & Projections - Left-join baseline and all three scenario-specific g_CO2e_1L columns into a single final_table. - Render as an interactive datatable or export to CSV.
Final Output (final_table): - Identifiers: Entity, Code, Region - Baseline Metrics: co2_intensity, g_CO2e_1L - 2030 Target Scenario: co2_intensity_2030, g_CO2e_1L_2030 - Half-Fossil Scenario: co2_intensity_half, g_CO2e_1L_half - Clean-Grid Scenario: co2_intensity_clean, g_CO2e_1L_clean
# 6.1 2030 target
grid_2030 <- grid_data_ci |>
select(-co2_intensity) |>
distinct(Code, .keep_all = TRUE) |>
left_join(
select(target_intensities, Code, co2_intensity_2030),
by = "Code"
) |>
rename(co2_intensity = co2_intensity_2030)
emissions_components_2030 <- grid_2030 |>
select(Entity = Code, co2_intensity) |>
crossing(systems_long) |>
mutate(
g_GWP_1L = if_else(
component == "prod",
g_GWP_1L * (co2_intensity / 28),
g_GWP_1L
)
) |>
group_by(Entity, system) |>
mutate(
percent_emissions = g_GWP_1L / sum(g_GWP_1L[component != "total"]) * 100
) |>
ungroup()
emissions_totals_2030 <- emissions_components_2030 |>
filter(component != "total") |>
group_by(Entity, system) |>
summarise(total_g_GWP_1L = sum(g_GWP_1L), .groups = "drop")
emissions_wide_2030 <- emissions_totals_2030 |>
mutate(system = paste0("GWP_", system)) |>
pivot_wider(names_from = system, values_from = total_g_GWP_1L)
country_emissions_2030 <- grid_2030 |>
select(Code) |>
left_join(emissions_wide_2030, by = c("Code" = "Entity")) |>
filter(!is.na(Code) & Code != "")
country_1l_2030 <- country_emissions_2030 |>
left_join(select(grid_2030, Code, Region), by = "Code") |>
left_join(oxygen_mix, by = c("Region" = "Entity")) |>
rowwise() |>
mutate(
all_gwp_na = if_all(all_of(system_cols), is.na),
g_CO2e_1L_2030 = if_else(
all_gwp_na,
NA_real_,
sum(
c_across(all_of(system_cols)) * c_across(all_of(share_cols)),
na.rm = TRUE
) /
100
)
) |>
ungroup() |>
select(Code, g_CO2e_1L_2030) |>
filter(!is.na(g_CO2e_1L_2030))
# 6.2 –50% fossil
grid_half <- grid_data_ci |>
select(-co2_intensity) |>
distinct(Code, .keep_all = TRUE) |>
left_join(
select(target_intensities, Code, co2_intensity_half),
by = "Code"
) |>
rename(co2_intensity = co2_intensity_half)
emissions_components_half <- grid_half |>
select(Entity = Code, co2_intensity) |>
crossing(systems_long) |>
mutate(
g_GWP_1L = if_else(
component == "prod",
g_GWP_1L * (co2_intensity / 28),
g_GWP_1L
)
) |>
group_by(Entity, system) |>
mutate(
percent_emissions = g_GWP_1L / sum(g_GWP_1L[component != "total"]) * 100
) |>
ungroup()
emissions_totals_half <- emissions_components_half |>
filter(component != "total") |>
group_by(Entity, system) |>
summarise(total_g_GWP_1L = sum(g_GWP_1L), .groups = "drop")
emissions_wide_half <- emissions_totals_half |>
mutate(system = paste0("GWP_", system)) |>
pivot_wider(names_from = system, values_from = total_g_GWP_1L)
country_emissions_half <- grid_half |>
select(Code) |>
left_join(emissions_wide_half, by = c("Code" = "Entity")) |>
filter(!is.na(Code) & Code != "")
country_1l_half <- country_emissions_half |>
left_join(select(grid_half, Code, Region), by = "Code") |>
left_join(oxygen_mix, by = c("Region" = "Entity")) |>
rowwise() |>
mutate(
all_gwp_na = if_all(all_of(system_cols), is.na),
g_CO2e_1L_half = if_else(
all_gwp_na,
NA_real_,
sum(
c_across(all_of(system_cols)) * c_across(all_of(share_cols)),
na.rm = TRUE
) /
100
)
) |>
ungroup() |>
select(Code, g_CO2e_1L_half) |>
filter(!is.na(g_CO2e_1L_half))
# 6.3 0% fossil
grid_clean <- grid_data_ci |>
select(-co2_intensity) |>
distinct(Code, .keep_all = TRUE) |>
left_join(
select(target_intensities, Code, co2_intensity_clean),
by = "Code"
) |>
rename(co2_intensity = co2_intensity_clean)
emissions_components_clean <- grid_clean |>
select(Entity = Code, co2_intensity) |>
crossing(systems_long) |>
mutate(
g_GWP_1L = if_else(
component == "prod",
g_GWP_1L * (co2_intensity / 28),
g_GWP_1L
)
) |>
group_by(Entity, system) |>
mutate(
percent_emissions = g_GWP_1L / sum(g_GWP_1L[component != "total"]) * 100
) |>
ungroup()
emissions_totals_clean <- emissions_components_clean |>
filter(component != "total") |>
group_by(Entity, system) |>
summarise(total_g_GWP_1L = sum(g_GWP_1L), .groups = "drop")
emissions_wide_clean <- emissions_totals_clean |>
mutate(system = paste0("GWP_", system)) |>
pivot_wider(names_from = system, values_from = total_g_GWP_1L)
country_emissions_clean <- grid_clean |>
select(Code) |>
left_join(emissions_wide_clean, by = c("Code" = "Entity")) |>
filter(!is.na(Code) & Code != "")
country_1l_clean <- country_emissions_clean |>
left_join(select(grid_clean, Code, Region), by = "Code") |>
left_join(oxygen_mix, by = c("Region" = "Entity")) |>
rowwise() |>
mutate(
all_gwp_na = if_all(all_of(system_cols), is.na),
g_CO2e_1L_clean = if_else(
all_gwp_na,
NA_real_,
sum(
c_across(all_of(system_cols)) * c_across(all_of(share_cols)),
na.rm = TRUE
) /
100
)
) |>
ungroup() |>
select(Code, g_CO2e_1L_clean) |>
filter(!is.na(g_CO2e_1L_clean))
# 6.4 Final assembly & export
final_table <- target_intensities |>
select(
Entity,
Code,
Region,
co2_intensity,
co2_intensity_2030,
co2_intensity_half,
co2_intensity_clean,
g_CO2e_1L
) |>
left_join(country_1l_2030, by = "Code") |>
left_join(country_1l_half, by = "Code") |>
left_join(country_1l_clean, by = "Code")
# ── 6.5 Render final_table as an interactive DataTable ────────────────────────
DT::datatable(
final_table,
extensions = 'Buttons',
options = list(
pageLength = 25,
scrollX = TRUE,
dom = 'Bfrtip',
buttons = c('csv', 'excel')
),
caption = '### Emissions per Liter Projections'
)